A Simple Model of Evolution with Variable System Size 
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A simple model of biological extinction with variable system size which exhibits a power-law 
distribution of extinction event sizes is presented. The model is a generalization of a model recently 
introduced by Newman (Proc. R. Soc. Lond. B263, 1605 (1996)). Both analytical and numerical 
analysis show that the exponent of the power-law distribution depends only marginally on the growth 
rate g at which new species enter the system and is equal to the one of the original model in the 
limit g oo. A critical growth rate gc, below which the system dies out, can be found. Under these 
model assumptions stable ecosystems can only exist if the regrowth of species is sufficiently fast. 

PACS numbers: 87.10. -l-e, 05.40. -Hj 

The fact that extinction events seem to be episodic 
on all scales, as noted by Raup has aroused much 
interest in the last few years. Throughout the history 
of life on earth there have been many small extinction 
events, but very big ones have happened only rarely. A 
histogram of the frequency of extinction events of differ- 
ent sizes indicates a power law distribution p(s) = s~'^, 
where s denotes the number of species that go extinct 
in one event and p(s) denotes the frequency of events of 
size s. 

There are two mechanisms to explain mass extinc- 
tions. On the one hand, it is argued that coevolution can 
drive large proportions of an ecosystem into extinction 
and produce extinction events on all scales. Ecosystems 
might drive themselves into a critical state in which a 
small change (e.g. the mutation of a single species) can 
trigger an "avalanche" that may span the whole system. 
For this kind of dynamic Bak et al. have coined the 
name Self-Organized Criticality. Several simple models 
of evolution exhibiting SOC have been proposed, among 
them models by Kauffman and Johnson Q|, Bak and 
Sneppen Q, Manrubia and Paczuski [||. 

On the other hand, it is argued that mass extinctions 
find their origin in external influences. That situation is 
modelled by some recent work of Newman . He used 
a model belonging to the new class of so called "coher- 
ent noise" models recently introduced by Newman and 
Sneppen 0. These models are clearly not SOC but they 
nevertheless show a power law distribution of avalanche 
sizes. Newman compared his model with the analysis of 
the fossil record performed by Raup. The exponent r 
close to 2 that arises in this model is in good agreement 
with the fossil record. Thus Newman came to the con- 
clusion that there is no evidence for SOC as the major 
driving force for extinction. 

It can be generally observed that the majority of the 
models for biological evolution and extinction up to now 
considered work with a fixed number of species. This is a 
major drawback since it is in clear contrast with the bio- 



logical reality. After a major extinction event, the num- 
ber of species in the ecosystem is significantly reduced, 
and the process of regrowth of new species can take a 
long time. The fossil record |^ shows that the process of 
growth of species is commonly interrupted by extinction 
events. 

To our knowledge, models with variable system size 
have only been studied by Vandewalle and Auslool |^ and 
by Head and Rodgers jl^ . But in both cases the models 
do not explain the distribution of extinction events seen 
in the fossil record. The model of Vandewalle and Aus- 
lool is a tree model that grows infinitely, while the model 
of Head and Rodgers reaches a steady-state in which no 
major extinctions occur. As far as we know, none of the 
models with variable system size up to now considered 
can explain the distribution of extinction events seen in 
the fossil record. 

But every mechanism proposed for the explanation of 
mass extinctions must 

• explain the distribution seen in the fossil record, 

• face the fact that the number of species is not con- 
stant, but is reduced significantly after a major ex- 
tinction event. 

A priori it is not at all clear if a mechanism producing 
a certain distribution of extinction events will show the 
same distribution when the constraint of a fixed system 
size is released. Therefore it is very important to study 
models with variable system size. 

We propose here a generalization to the coherent noise 
model used by Newman, where the refilling of the sys- 
tem is done in finite time. Newman's model is defined 
as follows. The system consists of N species, each pos- 
sessing a threshold Xi of tolerance against stress, cho- 
sen from a probablity distribution J3thrcsh(2^)- At each 
time-step, a stress 77 is generated at random with a dis- 
tribution Psticss('7), and all species with Xi < rj are re- 
moved from the system and immediately replaced with 
new ones. Furthermore, a small fraction / of the species 
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is chosen at random and given new thresholds. That cor- 
responds to a probability of / for every species to undergo 
spontaneous mutation. 

In our model the fraction of species with Xi < r] is re- 
moved permanently from the system, but in every time- 
step there is some growth of new species. 

Note that the generalized model, like the original one, 
does not include explicitly interaction between species. 
There are two reasons to justify this model assumption. 
Firstly, previous work has shown that the coherent 
noise dynamic is very strong and can dominate inter- 
action dynamic. Secondly, the investigation of a model 
without interaction, that can reproduce important fea- 
tures of the fossil record, helps to clarify the influence of 
species' interaction on mass extinctions. 

The amount of newly introduced species per time-step 
should be proportional to the number of already ex- 
isting species, with some constant of proportionality g 
(the growth rate). This gives an unbounded exponen- 
tial growth, which is in good agreement with the data of 
Benton However, since recourses on earth are finite, 
the growth of the species must be limited as well. There- 
fore, we believe it is justified to introduce a logistic factor 
(1 — ), where iVmax is the maximal number of species 
that can be sustained with the available resources. With 
this factor it is possible to work with a finite model. A 
few comments on the fact that in nature this A^max is 
probably not constant in time will be given later. 

For the above reasons we want our system to grow ac- 
cording to the differential equation 
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Since our model is discrete, in time as well as in the num- 
ber of species, instead of (|l|) we use the corresponding 
difference equation 
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where At is one simulation time-step (usually set equal 
to 1). As AA^ has to be an integer, we use the frac- 
tional part of AN as probability to round up or down. 
In the limit g ^ (which corresponds to At — > 0) Equa- 
tion (||) reduces to Equation (|l|). In the limit g ^ oo 
Equation (||) becomes A7V = A^max — N, which means 
that our model reduces to the original one in the limit of 
an infinite growth rate. 

Now we can formulate our model: We set At = 1. At 
every time step, a stress value rj is chosen and all species 
with Xi < T] are removed. Then, an amount AA^ of new 
species is introduced into the system. Finally, a fraction 
/ of the species is assigned new thresholds. 

A typical evolution of the system size N in time is pre- 
sented in Figure |l|. The process of growth of new species 
is constantly disrupted by small extinction events. From 



time to time, bigger events, which disturb the system sig- 
nificantly, occur. A plot of the distribution of extinction 
events (Figure ||) shows a power-law decrease. Variation 
of the growth rate over several orders of magnitude does 
change the exponent only slightly. 

We can explain the exponent of the power-law by ex- 
tending the analysis of Sneppen and Newman to our 
model. The probability that species leave a small inter- 
vall Ax of the time averaged distribution p{x) is propor- 
tional to (/ + Piaofc{x))p{x), where Pmove(;c) is the prob- 
ability that a species with threshold x is hit by stress. 
Let a be a variable that measures the "emptiness" of 
the system, i.e. a oc (1 — N/Nj^ax)- The rate at which 
the intervall dx is repopulated is then proportional to 
(/(I — a) + ga{l — a))pthTcshix) in the limit A< 0. 
In equilibrium the rates of species' loss and repopulation 
balance, and we find the master-equation 

if +Pmovc{x))pix) = (/(I - a) + .9(5(1 - a))pthrcsh(a:) . 

(3) 

Note that we had to replace a by its time-averaged value 
a and that we can always take the limit At — > in the 
steady-state. After rearranging Equation (^ , we find 
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Equation (Q) can be solved if we choose how to normalize 
p{x) and a. Since we can think of the system as contain- 
ing A'ljiax species at any time step, from which there are 
N active and Nj^ax — N dead, it makes sense to normalize 
the sum of a and p{x) to unity, viz. 



1 = a -I- J p{x)dx . 



(5) 



That implies, nevertheless, that we do not normalize p{x) 
to unity. Rather, J p(x)dx gives the ratio N/N„^ax- 

For a we find, apart from the trivial solution a — 1, 
the solution a = {A — f)/g, with 
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For p(x), we find 



p{x) =A[1 
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We thus have the interesting result that apart from the 
overall factor 1 — a, which determines the average sys- 
tem size, the shape of p{x) is identical to that found by 
Sneppen and Newman. Since only the shape p{x), but 
not the overall factor, is responsible for the power-law 
distribution of extinction events (for details see 0) we 
find that, within the time averaged approximation, the 
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exponent r of the power-law decrease is exactly the same 
as in the original model, even for very small g. 

If we take the limit 5 oo in Equation (|^) we can 
restore the expression found by Sneppen and Newman, 
which was to be expected since our model reduces to the 
original one in that limit. In the region of very small g, 
we can read off from Equation that the system breaks 
down at a critical growth rate gc — A — f. This is the 
case when the growth rate is so small that the regrowth 
of species cannot compensate the successive extinction 
events. Every system with g < gc will eventually end up 
with = 0, regardless of the number of species at the 
beginning of the simulation. 

For the simulation results presented here we have used 
exponentially distributed stress only, i.e., Pstrcssl??) = 
exp{—r]/a)/<7. We did simulations with A^max between 
1000 and 10000. Figure || shows the dependence of the 
average system size N of g. We can clearly see the break- 
down of the system at gc- A measurement of the timc- 
averaged distribution of thresholds p{x) is presented in 
Figure ^. The exponent r of the power-law distribution 
of extinction events is found to be r = 1.9±0.1 for g = 10, 
T = 2.0±0.1 for g = 0.002, r = 2.05±0.1 for g = 4 x 10"^ 
(for exponentially distributed stress, a = 0.05, / = 10~^, 
Figure ||) . The exponent decreases slightly with increas- 
ing g. For g = 10, we have already good agreement 
with the exponent found by Newman and Sneppen |^] 
for g = 00, viz. t = 1.85 ± 0.03. 

An interesting feature of the original model by New- 
man and Sneppen is the existence of aftershocks, a series 
of smaller events following a large one. These aftershocks 
have their origin in the fact that after a large event the in- 
troduction of new species reduces significantly the mean 
threshold value, thus increasing the probability to get fur- 
ther events. Since the existence of aftershocks is a result 
of the immediate refilling of the system after an event, we 
cannot necessarily expect to see aftershocks when the re- 
filling is done in finite time, especially at a small growth 
rate. Numerical simulations show that there are after- 
shocks for larger values of g, but when g approaches gc, 
aftershocks cannot clearly be identified anymore. The 
region where this happens is that in which the average 
system size decreases rapidly with g. For these values of 
g, the typical time the system needs to regrow the amount 
of species lost in a major event exceeds the typical time 
needed to create a major stress value. In Figure H the 
region in which we do not find aftershocks is between 
g ~ gc = 1.3 X 10^^ and about g = 5 x 10^*. A typical 
example for a series of events in a system with g close to 
gc is presented in Figure ^. 

Sneppen and Newman argued that the existence of 
aftershocks might provide a measure to distinguish be- 
tween coherent-noise driven systems and SOC systems. 
This is certainly true in the sense that systems exhibit- 
ing aftershocks are better candidates for coherent-noise 
driven systems rather than for SOC systems. But our 



simulations show that there are systems without clear af- 
tershocks that still should be classified as coherent-noise 
driven. 

We have focused on logistic growth since we believe it 
is suitable for the study of mass extinctions. In principle 
it is possible to use different types of growth. We have 
done some simulations with linear growth, where in ev- 
ery time-step a fixed amount of new species is introduced 
into the system, as long as < A^max- These simulations 
indicate that the respective type of growth used does not 
affect the appearance of a power-law distribution with 
exponent almost independent from the growth rate. But 
whether aftershocks appear or not, is indeed dependend 
on the type of growth we choose. In a system with lin- 
ear growth aftershocks can be seen clearly even for small 
growth rates. 

If we want to use a coherent noise model with variable 
system size as a model of biological evolution, some re- 
marks about the meaning of A^max are necessary. The 
fact of allowing the regrowth of species in finite time, 
instead of refilling the system immediately, represents a 
first step closer to reality. But for ecosystems it is cer- 
tainly not a good assumption to keep the maximal system 
size A'max fixed, since the number of species an ecosys- 
tem can contain depends on the interaction of species 
themselves. Therefore, a next step could be to change 
A'max after every extinction, e.g., up or down by chance 
and by an amount proportional to the size of the event. 
This is motivated by the fact that bigger events are ex- 
pected to be correlated with a more profound restruc- 
turing of the ecosystem, and as simulations show we still 
find power-law distributions with exponents r w 2. The 
behaviour of such a system has a very rich structure with 
long times of relatively little change (stasis) and sudden 
bursts of evolutionary activity (punctuated equilibrium), 
where a major extinction event is followed by a regrowth 
of species to a system size much bigger than the one be- 
fore the event. The so found curves of the system size A^ 
agree qualitatively well with the fossil record 

We have generalized a coherent noise model to a model 
with variable system size. The most important feature 
of coherent noise models, the power-law distribution of 
event sizes with an exponent close to 2, does not change 
under the generalization. This means that the validity of 
Newman's approach to explain biological extinction with 
a coherent noise model is not affected by the regrowth 
of species in finite time. An interesting new feature that 
emerges from a variable system size is the existence of 
a critical growth rate gc- Systems with g < gc will al- 
ways end up with A^ = after some time. Therefore in 
a world in which the regrowth of species is too slow to 
compensate external influences no stable ecosystems can 
exist. In the framework of our model we conclude that 
the process of mutation and diversification of species at 
sufficiently high rate is necessary for the stability of life 
on earth. 
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extinction event size s 
FIG. 2. The distribution of extinction events for a sys- 
tem with exponentially distributed stress, a = 0.05 and 
Amax = 10000. The growth rate is, from bottom to top, 
3 = 4 X 10~^, g = 0.002, g = 10. It can be seen that 
the power-law behavior does depend only marginally on the 
growth rate. The curves have been rescaled so as not to over- 
lap. 
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FIG. 3. The average system size iV versus the growth 
rate g. We used exponentially distributed stress with a — 0.05 
and / — 10~^. The solid line is the analytic expression, the 
points are the simulation results. 
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FIG. 1. The evolution of the system size A in time. 

parameters are <; = 4 x 10~^, a = 0.05, / = 10~^, 
Amax = 1000 with exponentially distributed stress. 
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FIG. 4. The time averaged distribution p{x). The param- 
eters used are g = 0.002, a = 0.05, and / = 5 x 10~* with 
exponentially distributed stress. The solid line is the analytic 
expression, the points are the simulation results. 
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FIG. 5. A series of extinction events. The parameters used 
are g = 4 X 10~®, a = 0.05, and / = 5 x 10""* with ex- 
ponentially distributed stress. Aftershocks cannot clearly be 
identified. 
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